library(plyr)
library(igraph)
library(RgoogleMaps)
library(RColorBrewer)
library(fmsb)
library(Hmisc)
library(MASS)
library(stargazer)
library(tidyverse)
library(multcomp)
library(WDI)
library(countrycode)

setwd("~/Dropbox/China Huawei Paper/Production Materials/Replication/")
master.data <- read.csv("dataFile_China_Huawei.csv"); head(master.data)
savedir <- '~/Dropbox/China Huawei Paper/Production Materials/Figures/'

# V-Dem addition
vdem <- read.csv("V-Dem-CY-Full+Others-v12.csv"); head(vdem); dim(vdem)
vdem <- vdem[,c(which(colnames(vdem)=="COWcode"), which(colnames(vdem)=="year"), 
                which(colnames(vdem)=="v2xcl_rol"), which(colnames(vdem)=="v2xlg_legcon"), which(colnames(vdem)=="v2xel_frefair"), which(colnames(vdem)=="v2meharjrn")
                )]; head(vdem)
vdem <- vdem[vdem$year>=2000,]
vdem$cabb <- countrycode(sourcevar=vdem$COWcode, origin="cown", destination="cowc")
vdem$COWcode=NULL
nrow(master.data)
master.data <- merge(master.data, vdem, by=c("cabb","year"), all.x=T); nrow(master.data)


# 1,000,000
dat <- master.data
dat$treated = 0
threshold=1000000

dat = as.data.frame(dat %>%
                      group_by(cabb) %>%
                      dplyr::mutate(treated = ifelse(row_number() == 1, 0, NA),
                                    treated = ifelse(china.huawei > threshold, 1, treated)) %>%
                      fill(treated))
dat$treated[is.na(dat$treated)==T] <- 0
dat = dat[dat$v2x_polyarchy>0.42,]

dat$v2x_frassoc_thick_treated = dat$treated*dat$v2x_frassoc_thick
dat$v2xcl_rol_treated = dat$treated*dat$v2xcl_rol
dat$v2xlg_legcon_treated = dat$treated*dat$v2xlg_legcon
dat$v2xel_frefair_treated = dat$treated*dat$v2xel_frefair
dat$v2meharjrn_treated = dat$treated*dat$v2meharjrn


# Lagging variables
paneldata.lags <- function(X, unitvar, timevar, lagvars, lags=c(1:3)){
  stopifnot(length(lagvars)>=1)
  X <- X[order(X[,timevar]),] # just in case it's not sorted.
  
  innertask <- function(Y, lagvars, lags) {     
    E <- labels <- NULL
    for (v in lagvars) {
      for (i in lags) {E <- cbind(E, Lag(Y[,v], i))}
      labels <- c(labels, paste(v, ".l", lags, sep=""))
    }
    colnames(E) <- labels
    cbind(Y, E)
  }
  do.call("rbind", by(X, X[,unitvar], innertask, lagvars, lags)) 
}
dat <- dat[order(dat$cabb, dat$year, decreasing=F),]
dat <- paneldata.lags(X=dat, unitvar="cabb", timevar="year", lagvars=c("treated", "v2x_frassoc_thick", "v2xel_frefair", "v2xcl_rol", "v2xlg_legcon", "v2meharjrn"))


####### Freedom of Association

out1 = lm(v2smgovfilprc ~ treated + v2x_frassoc_thick_treated + treated.l1 + treated.l2 + treated.l3 
          + v2x_frassoc_thick + v2x_frassoc_thick.l1 + v2x_frassoc_thick.l2 + v2x_frassoc_thick.l3
          + pres.election + coup.attempts + successful.coups + icews.protest + icews.repression + GDPcur + pcGDPcur + electricitypc + as.factor(cabb), data=dat)
summary(out1)

out2 = lm(v2smgovshut ~ treated + v2x_frassoc_thick_treated + treated.l1 + treated.l2 + treated.l3 
          + v2x_frassoc_thick + v2x_frassoc_thick.l1 + v2x_frassoc_thick.l2 + v2x_frassoc_thick.l3
          + pres.election + coup.attempts + successful.coups + icews.protest + icews.repression + GDPcur + pcGDPcur + electricitypc + as.factor(cabb), data=dat)
summary(out2)

out3 = lm(v2smgovsmmon ~ treated + v2x_frassoc_thick_treated + treated.l1 + treated.l2 + treated.l3 
          + v2x_frassoc_thick + v2x_frassoc_thick.l1 + v2x_frassoc_thick.l2 + v2x_frassoc_thick.l3
          + pres.election + coup.attempts + successful.coups + icews.protest + icews.repression + GDPcur + pcGDPcur + electricitypc + as.factor(cabb), data=dat)
summary(out3)

out4 = lm(v2smarrest ~ treated + v2x_frassoc_thick_treated + treated.l1 + treated.l2 + treated.l3 
          + v2x_frassoc_thick + v2x_frassoc_thick.l1 + v2x_frassoc_thick.l2 + v2x_frassoc_thick.l3
          + pres.election + coup.attempts + successful.coups + icews.protest + icews.repression + GDPcur + pcGDPcur + electricitypc + as.factor(cabb), data=dat)
summary(out4)

stargazer(out1,out2,out3,out4,
          add.lines=list(c("Controls"  ,"$\\checkmark$","$\\checkmark$","$\\checkmark$","$\\checkmark$"),
                         c("Treatment Lags"  ,"$\\checkmark$","$\\checkmark$","$\\checkmark$","$\\checkmark$"),
                         c("Subcomponent Lags"  ,"$\\checkmark$","$\\checkmark$","$\\checkmark$","$\\checkmark$"),
                         c("Country Fixed Effects"  ,"$\\checkmark$","$\\checkmark$","$\\checkmark$","$\\checkmark$")
          ),
          title="How Democratic Institutions Moderate the Effect of Huawei Transfers",
          label="Table: Vdem Subcomponents",
          omit=c("cabb",".l1",".l2",".l3","pres.election","coup.attempts","successful.coups","icews.protest","icews.repression","GDPcur","electricitypc"),
          font.size="footnotesize",
          dep.var.labels = c("Internet Filtering","Internet Shutdowns","Social Media Monitoring","Arrests for Political Content"),
          covariate.labels = c("Treatment", "Freedom of Association $\\times$ Treatment", "Freedom of Association"
          )
)

frassoc.fil <- rbind(marginal.effects(out1, conf=95, beta.coef=2, interaction=T, interaction.coef=3, Z=min(dat$v2x_frassoc_thick))[2,], marginal.effects(out1, conf=95, beta.coef=2, interaction=T, interaction.coef=3, Z=max(dat$v2x_frassoc_thick))[2,])

frassoc.ish <- rbind(marginal.effects(out2, conf=95, beta.coef=2, interaction=T, interaction.coef=3, Z=min(dat$v2x_frassoc_thick))[2,], marginal.effects(out2, conf=95, beta.coef=2, interaction=T, interaction.coef=3, Z=max(dat$v2x_frassoc_thick))[2,])

frassoc.smm <- rbind(marginal.effects(out3, conf=95, beta.coef=2, interaction=T, interaction.coef=3, Z=min(dat$v2x_frassoc_thick))[2,], marginal.effects(out3, conf=95, beta.coef=2, interaction=T, interaction.coef=3, Z=max(dat$v2x_frassoc_thick))[2,])

frassoc.arr <- rbind(marginal.effects(out4, conf=95, beta.coef=2, interaction=T, interaction.coef=3, Z=min(dat$v2x_frassoc_thick))[2,], marginal.effects(out4, conf=95, beta.coef=2, interaction=T, interaction.coef=3, Z=max(dat$v2x_frassoc_thick))[2,])


####### Rule of Law

out1 = lm(v2smgovfilprc ~ treated + v2xcl_rol_treated + treated.l1 + treated.l2 + treated.l3 
          + v2xcl_rol + v2xcl_rol.l1 + v2xcl_rol.l2 + v2xcl_rol.l3
          + pres.election + coup.attempts + successful.coups + icews.protest + icews.repression + GDPcur + pcGDPcur + electricitypc + as.factor(cabb), data=dat)
summary(out1)

out2 = lm(v2smgovshut ~ treated + v2xcl_rol_treated + treated.l1 + treated.l2 + treated.l3 
          + v2xcl_rol + v2xcl_rol.l1 + v2xcl_rol.l2 + v2xcl_rol.l3
          + pres.election + coup.attempts + successful.coups + icews.protest + icews.repression + GDPcur + pcGDPcur + electricitypc + as.factor(cabb), data=dat)
summary(out2)

out3 = lm(v2smgovsmmon ~ treated + v2xcl_rol_treated + treated.l1 + treated.l2 + treated.l3 
          + v2xcl_rol + v2xcl_rol.l1 + v2xcl_rol.l2 + v2xcl_rol.l3
          + pres.election + coup.attempts + successful.coups + icews.protest + icews.repression + GDPcur + pcGDPcur + electricitypc + as.factor(cabb), data=dat)
summary(out3)

out4 = lm(v2smarrest ~ treated + v2xcl_rol_treated + treated.l1 + treated.l2 + treated.l3 
          + v2xcl_rol + v2xcl_rol.l1 + v2xcl_rol.l2 + v2xcl_rol.l3
          + pres.election + coup.attempts + successful.coups + icews.protest + icews.repression + GDPcur + pcGDPcur + electricitypc + as.factor(cabb), data=dat)
summary(out4)

stargazer(out1,out2,out3,out4,
          add.lines=list(c("Controls"  ,"$\\checkmark$","$\\checkmark$","$\\checkmark$","$\\checkmark$"),
                         c("Treatment Lags"  ,"$\\checkmark$","$\\checkmark$","$\\checkmark$","$\\checkmark$"),
                         c("Subcomponent Lags"  ,"$\\checkmark$","$\\checkmark$","$\\checkmark$","$\\checkmark$"),
                         c("Country Fixed Effects"  ,"$\\checkmark$","$\\checkmark$","$\\checkmark$","$\\checkmark$")
          ),
          title="How Democratic Institutions Moderate the Effect of Huawei Transfers",
          label="Table: Vdem Subcomponents",
          omit=c("cabb",".l1",".l2",".l3","pres.election","coup.attempts","successful.coups","icews.protest","icews.repression","GDPcur","electricitypc"),
          font.size="footnotesize",
          dep.var.labels = c("Internet Filtering","Internet Shutdowns","Social Media Monitoring","Arrests for Political Content"), 
          covariate.labels = c("Treatment", "Rule of Law $\\times$ Treatment", "Rule of Law"
          )
)

rol.fil <- rbind(marginal.effects(out1, conf=95, beta.coef=2, interaction=T, interaction.coef=3, Z=min(dat$v2xcl_rol))[2,], marginal.effects(out1, conf=95, beta.coef=2, interaction=T, interaction.coef=3, Z=max(dat$v2xcl_rol))[2,])

rol.ish <- rbind(marginal.effects(out2, conf=95, beta.coef=2, interaction=T, interaction.coef=3, Z=min(dat$v2xcl_rol))[2,], marginal.effects(out2, conf=95, beta.coef=2, interaction=T, interaction.coef=3, Z=max(dat$v2xcl_rol))[2,])

rol.smm <- rbind(marginal.effects(out3, conf=95, beta.coef=2, interaction=T, interaction.coef=3, Z=min(dat$v2xcl_rol))[2,], marginal.effects(out3, conf=95, beta.coef=2, interaction=T, interaction.coef=3, Z=max(dat$v2xcl_rol))[2,])

rol.arr <- rbind(marginal.effects(out4, conf=95, beta.coef=2, interaction=T, interaction.coef=3, Z=min(dat$v2xcl_rol))[2,], marginal.effects(out4, conf=95, beta.coef=2, interaction=T, interaction.coef=3, Z=max(dat$v2xcl_rol))[2,])


####### Legislative Constraints

out1 = lm(v2smgovfilprc ~ treated + v2xlg_legcon_treated + treated.l1 + treated.l2 + treated.l3 
          + v2xlg_legcon + v2xlg_legcon.l1 + v2xlg_legcon.l2 + v2xlg_legcon.l3
          + pres.election + coup.attempts + successful.coups + icews.protest + icews.repression + GDPcur + pcGDPcur + electricitypc + as.factor(cabb), data=dat)
summary(out1)

out2 = lm(v2smgovshut ~ treated + v2xlg_legcon_treated + treated.l1 + treated.l2 + treated.l3 
          + v2xlg_legcon + v2xlg_legcon.l1 + v2xlg_legcon.l2 + v2xlg_legcon.l3
          + pres.election + coup.attempts + successful.coups + icews.protest + icews.repression + GDPcur + pcGDPcur + electricitypc + as.factor(cabb), data=dat)
summary(out2)

out3 = lm(v2smgovsmmon ~ treated + v2xlg_legcon_treated + treated.l1 + treated.l2 + treated.l3 
          + v2xlg_legcon + v2xlg_legcon.l1 + v2xlg_legcon.l2 + v2xlg_legcon.l3
          + pres.election + coup.attempts + successful.coups + icews.protest + icews.repression + GDPcur + pcGDPcur + electricitypc + as.factor(cabb), data=dat)
summary(out3)

out4 = lm(v2smarrest ~ treated + v2xlg_legcon_treated + treated.l1 + treated.l2 + treated.l3 
          + v2xlg_legcon + v2xlg_legcon.l1 + v2xlg_legcon.l2 + v2xlg_legcon.l3
          + pres.election + coup.attempts + successful.coups + icews.protest + icews.repression + GDPcur + pcGDPcur + electricitypc + as.factor(cabb), data=dat)
summary(out4)

stargazer(out1,out2,out3,out4,
          add.lines=list(c("Controls"  ,"$\\checkmark$","$\\checkmark$","$\\checkmark$","$\\checkmark$"),
                         c("Treatment Lags"  ,"$\\checkmark$","$\\checkmark$","$\\checkmark$","$\\checkmark$"),
                         c("Subcomponent Lags"  ,"$\\checkmark$","$\\checkmark$","$\\checkmark$","$\\checkmark$"),
                         c("Country Fixed Effects"  ,"$\\checkmark$","$\\checkmark$","$\\checkmark$","$\\checkmark$")
          ),
          title="How Democratic Institutions Moderate the Effect of Huawei Transfers",
          label="Table: Vdem Subcomponents",
          omit=c("cabb",".l1",".l2",".l3","pres.election","coup.attempts","successful.coups","icews.protest","icews.repression","GDPcur","electricitypc"),
          font.size="footnotesize",
          dep.var.labels = c("Internet Filtering","Internet Shutdowns","Social Media Monitoring","Arrests for Political Content"), 
          covariate.labels = c("Treatment", "Legislative Constraints $\\times$ Treatment", "Legislative Constraints"
          )
)

legcon.fil <- rbind(marginal.effects(out1, conf=95, beta.coef=2, interaction=T, interaction.coef=3, Z=min(dat$v2xlg_legcon))[2,], marginal.effects(out1, conf=95, beta.coef=2, interaction=T, interaction.coef=3, Z=max(dat$v2xlg_legcon))[2,])

legcon.ish <- rbind(marginal.effects(out2, conf=95, beta.coef=2, interaction=T, interaction.coef=3, Z=min(dat$v2xlg_legcon))[2,], marginal.effects(out2, conf=95, beta.coef=2, interaction=T, interaction.coef=3, Z=max(dat$v2xlg_legcon))[2,])

legcon.smm <- rbind(marginal.effects(out3, conf=95, beta.coef=2, interaction=T, interaction.coef=3, Z=min(dat$v2xlg_legcon))[2,], marginal.effects(out3, conf=95, beta.coef=2, interaction=T, interaction.coef=3, Z=max(dat$v2xlg_legcon))[2,])

legcon.arr <- rbind(marginal.effects(out4, conf=95, beta.coef=2, interaction=T, interaction.coef=3, Z=min(dat$v2xlg_legcon))[2,], marginal.effects(out4, conf=95, beta.coef=2, interaction=T, interaction.coef=3, Z=max(dat$v2xlg_legcon))[2,])


####### Free and Fair Elections

out1 = lm(v2smgovfilprc ~ treated + v2xel_frefair_treated + treated.l1 + treated.l2 + treated.l3 
          + v2xel_frefair + v2xel_frefair.l1 + v2xel_frefair.l2 + v2xel_frefair.l3
          + pres.election + coup.attempts + successful.coups + icews.protest + icews.repression + GDPcur + pcGDPcur + electricitypc + as.factor(cabb), data=dat)
summary(out1)

out2 = lm(v2smgovshut ~ treated + v2xel_frefair_treated + treated.l1 + treated.l2 + treated.l3 
          + v2xel_frefair + v2xel_frefair.l1 + v2xel_frefair.l2 + v2xel_frefair.l3
          + pres.election + coup.attempts + successful.coups + icews.protest + icews.repression + GDPcur + pcGDPcur + electricitypc + as.factor(cabb), data=dat)
summary(out2)

out3 = lm(v2smgovsmmon ~ treated + v2xel_frefair_treated + treated.l1 + treated.l2 + treated.l3 
          + v2xel_frefair + v2xel_frefair.l1 + v2xel_frefair.l2 + v2xel_frefair.l3
          + pres.election + coup.attempts + successful.coups + icews.protest + icews.repression + GDPcur + pcGDPcur + electricitypc + as.factor(cabb), data=dat)
summary(out3)

out4 = lm(v2smarrest ~ treated + v2xel_frefair_treated + treated.l1 + treated.l2 + treated.l3 
          + v2xel_frefair + v2xel_frefair.l1 + v2xel_frefair.l2 + v2xel_frefair.l3
          + pres.election + coup.attempts + successful.coups + icews.protest + icews.repression + GDPcur + pcGDPcur + electricitypc + as.factor(cabb), data=dat)
summary(out4)

stargazer(out1,out2,out3,out4,
          add.lines=list(c("Controls"  ,"$\\checkmark$","$\\checkmark$","$\\checkmark$","$\\checkmark$"),
                         c("Treatment Lags"  ,"$\\checkmark$","$\\checkmark$","$\\checkmark$","$\\checkmark$"),
                         c("Subcomponent Lags"  ,"$\\checkmark$","$\\checkmark$","$\\checkmark$","$\\checkmark$"),
                         c("Country Fixed Effects"  ,"$\\checkmark$","$\\checkmark$","$\\checkmark$","$\\checkmark$")
          ),
          title="How Democratic Institutions Moderate the Effect of Huawei Transfers",
          label="Table: Vdem Subcomponents",
          omit=c("cabb",".l1",".l2",".l3","pres.election","coup.attempts","successful.coups","icews.protest","icews.repression","GDPcur","electricitypc"),
          font.size="footnotesize",
          dep.var.labels = c("Internet Filtering","Internet Shutdowns","Social Media Monitoring","Arrests for Political Content"), 
          covariate.labels = c("Treatment", "Free and Fair Elections $\\times$ Treatment", "Free and Fair Elections"
          )
)

frefair.fil <- rbind(marginal.effects(out1, conf=95, beta.coef=2, interaction=T, interaction.coef=3, Z=min(dat$v2xel_frefair))[2,], marginal.effects(out1, conf=95, beta.coef=2, interaction=T, interaction.coef=3, Z=max(dat$v2xel_frefair))[2,])

frefair.ish <- rbind(marginal.effects(out2, conf=95, beta.coef=2, interaction=T, interaction.coef=3, Z=min(dat$v2xel_frefair))[2,], marginal.effects(out2, conf=95, beta.coef=2, interaction=T, interaction.coef=3, Z=max(dat$v2xel_frefair))[2,])

frefair.smm <- rbind(marginal.effects(out3, conf=95, beta.coef=2, interaction=T, interaction.coef=3, Z=min(dat$v2xel_frefair))[2,], marginal.effects(out3, conf=95, beta.coef=2, interaction=T, interaction.coef=3, Z=max(dat$v2xel_frefair))[2,])

frefair.arr <- rbind(marginal.effects(out4, conf=95, beta.coef=2, interaction=T, interaction.coef=3, Z=min(dat$v2xel_frefair))[2,], marginal.effects(out4, conf=95, beta.coef=2, interaction=T, interaction.coef=3, Z=max(dat$v2xel_frefair))[2,])


####### Journalist harassment

out1 = lm(v2smgovfilprc ~ treated + v2meharjrn_treated + treated.l1 + treated.l2 + treated.l3 
          + v2meharjrn + v2meharjrn.l1 + v2meharjrn.l2 + v2meharjrn.l3
          + pres.election + coup.attempts + successful.coups + icews.protest + icews.repression + GDPcur + pcGDPcur + electricitypc + as.factor(cabb), data=dat)
summary(out1)

out2 = lm(v2smgovshut ~ treated + v2meharjrn_treated + treated.l1 + treated.l2 + treated.l3 
          + v2meharjrn + v2meharjrn.l1 + v2meharjrn.l2 + v2meharjrn.l3
          + pres.election + coup.attempts + successful.coups + icews.protest + icews.repression + GDPcur + pcGDPcur + electricitypc + as.factor(cabb), data=dat)
summary(out2)

out3 = lm(v2smgovsmmon ~ treated + v2meharjrn_treated + treated.l1 + treated.l2 + treated.l3 
          + v2meharjrn + v2meharjrn.l1 + v2meharjrn.l2 + v2meharjrn.l3
          + pres.election + coup.attempts + successful.coups + icews.protest + icews.repression + GDPcur + pcGDPcur + electricitypc + as.factor(cabb), data=dat)
summary(out3)

out4 = lm(v2smarrest ~ treated + v2meharjrn_treated + treated.l1 + treated.l2 + treated.l3 
          + v2meharjrn + v2meharjrn.l1 + v2meharjrn.l2 + v2meharjrn.l3
          + pres.election + coup.attempts + successful.coups + icews.protest + icews.repression + GDPcur + pcGDPcur + electricitypc + as.factor(cabb), data=dat)
summary(out4)

stargazer(out1,out2,out3,out4,
          add.lines=list(c("Controls"  ,"$\\checkmark$","$\\checkmark$","$\\checkmark$","$\\checkmark$"),
                         c("Treatment Lags"  ,"$\\checkmark$","$\\checkmark$","$\\checkmark$","$\\checkmark$"),
                         c("Subcomponent Lags"  ,"$\\checkmark$","$\\checkmark$","$\\checkmark$","$\\checkmark$"),
                         c("Country Fixed Effects"  ,"$\\checkmark$","$\\checkmark$","$\\checkmark$","$\\checkmark$")
          ),
          title="How Democratic Institutions Moderate the Effect of Huawei Transfers",
          label="Table: Vdem Subcomponents",
          omit=c("cabb",".l1",".l2",".l3","pres.election","coup.attempts","successful.coups","icews.protest","icews.repression","GDPcur","electricitypc"),
          font.size="footnotesize",
          dep.var.labels = c("Internet Filtering","Internet Shutdowns","Social Media Monitoring","Arrests for Political Content"), 
          covariate.labels = c("Treatment", "Journalist Harassment $\\times$ Treatment", "Journalist Harassment"
          )
)

meharjrn.fil <- rbind(marginal.effects(out1, conf=95, beta.coef=2, interaction=T, interaction.coef=3, Z=min(dat$v2meharjrn))[2,], marginal.effects(out1, conf=95, beta.coef=2, interaction=T, interaction.coef=3, Z=max(dat$v2meharjrn))[2,])

meharjrn.ish <- rbind(marginal.effects(out2, conf=95, beta.coef=2, interaction=T, interaction.coef=3, Z=min(dat$v2meharjrn))[2,], marginal.effects(out2, conf=95, beta.coef=2, interaction=T, interaction.coef=3, Z=max(dat$v2meharjrn))[2,])

meharjrn.smm <- rbind(marginal.effects(out3, conf=95, beta.coef=2, interaction=T, interaction.coef=3, Z=min(dat$v2meharjrn))[2,], marginal.effects(out3, conf=95, beta.coef=2, interaction=T, interaction.coef=3, Z=max(dat$v2meharjrn))[2,])

meharjrn.arr <- rbind(marginal.effects(out4, conf=95, beta.coef=2, interaction=T, interaction.coef=3, Z=min(dat$v2meharjrn))[2,], marginal.effects(out4, conf=95, beta.coef=2, interaction=T, interaction.coef=3, Z=max(dat$v2meharjrn))[2,])




# Internet Filtering
models=5
adjust=0.1
par(mar=c(bottom=4.0,left=6.5,top=1.0,right=0.2), lheight=0.8,xpd=F) # xpd lets plot outside graph area
plot(-1,-1, type="p", axes=F, xlab="", ylab="", xlim=c(-2,2), ylim=c(0.5,to=models+0.5), cex.lab=1.3, cex.main=1.3, main="Internet Filtering")
axis(1, at=seq(-2,2,by=0.5), line=0, las=0, tick=T, cex.axis=1.3)
axis(2, at=seq(1,models, by=1), label=rev(c("Freedom of\nAssociation", "Rule of Law", "Legislative\nConstraints", "Free & Fair\nElections", "Journalist\nHarassment")), line=0, col="Black", col.axis="Black", col.ticks="Black", cex.axis=1.2, tick=F, las=2)
abline(v=0, lty=2)
legend("topright",bty="n",lwd=3,col=c("darkgray","black"),c("Maximum Observed Value","Minimum Observed Value"), cex=1.1)

unit=1
segments(x0=meharjrn.fil[1,2], y0=unit-adjust, x1=meharjrn.fil[1,3], y1=unit-adjust, lwd=1.5, bg="black", col="black")
points(x=meharjrn.fil[1,1], y=unit-adjust, pch=21, bg="black", col="black", cex=1.5)
segments(x0=meharjrn.fil[2,2], y0=unit+adjust, x1=meharjrn.fil[2,3], y1=unit+adjust, lwd=1.5, bg="darkgray", col="darkgray")
points(x=meharjrn.fil[2,1], y=unit+adjust, pch=21, bg="darkgray", col="darkgray", cex=1.5)

unit=2
segments(x0=frefair.fil[1,2], y0=unit-adjust, x1=frefair.fil[1,3], y1=unit-adjust, lwd=1.5, bg="black", col="black")
points(x=frefair.fil[1,1], y=unit-adjust, pch=21, bg="black", col="black", cex=1.5)
segments(x0=frefair.fil[2,2], y0=unit+adjust, x1=frefair.fil[2,3], y1=unit+adjust, lwd=1.5, bg="darkgray", col="darkgray")
points(x=frefair.fil[2,1], y=unit+adjust, pch=21, bg="darkgray", col="darkgray", cex=1.5)

unit=3
segments(x0=legcon.fil[1,2], y0=unit-adjust, x1=legcon.fil[1,3], y1=unit-adjust, lwd=1.5, bg="black", col="black")
points(x=legcon.fil[1,1], y=unit-adjust, pch=21, bg="black", col="black", cex=1.5)
segments(x0=legcon.fil[2,2], y0=unit+adjust, x1=legcon.fil[2,3], y1=unit+adjust, lwd=1.5, bg="darkgray", col="darkgray")
points(x=legcon.fil[2,1], y=unit+adjust, pch=21, bg="darkgray", col="darkgray", cex=1.5)

unit=4
segments(x0=rol.fil[1,2], y0=unit-adjust, x1=rol.fil[1,3], y1=unit-adjust, lwd=1.5, bg="black", col="black")
points(x=rol.fil[1,1], y=unit-adjust, pch=21, bg="black", col="black", cex=1.5)
segments(x0=rol.fil[2,2], y0=unit+adjust, x1=rol.fil[2,3], y1=unit+adjust, lwd=1.5, bg="darkgray", col="darkgray")
points(x=rol.fil[2,1], y=unit+adjust, pch=21, bg="darkgray", col="darkgray", cex=1.5)

unit=5
segments(x0=frassoc.fil[1,2], y0=unit-adjust, x1=frassoc.fil[1,3], y1=unit-adjust, lwd=1.5, bg="black", col="black")
points(x=frassoc.fil[1,1], y=unit-adjust, pch=21, bg="black", col="black", cex=1.5)
segments(x0=frassoc.fil[2,2], y0=unit+adjust, x1=frassoc.fil[2,3], y1=unit+adjust, lwd=1.5, bg="darkgray", col="darkgray")
points(x=frassoc.fil[2,1], y=unit+adjust, pch=21, bg="darkgray", col="darkgray", cex=1.5)

dev.print(file=paste(savedir,'backsliding_Filtering.pdf',sep=''), device=pdf)

# Internet Shutdowns
models=5
adjust=0.1
par(mar=c(bottom=4.0,left=6.5,top=1.0,right=0.2), lheight=0.8,xpd=F) # xpd lets plot outside graph area
plot(-1,-1, type="p", axes=F, xlab="", ylab="", xlim=c(-2,2), ylim=c(0.5,to=models+0.5), cex.lab=1.3, cex.main=1.3, main="Internet Shutdowns")
axis(1, at=seq(-2,2,by=0.5), line=0, las=0, tick=T, cex.axis=1.3)
axis(2, at=seq(1,models, by=1), label=rev(c("Freedom of\nAssociation", "Rule of Law", "Legislative\nConstraints", "Free & Fair\nElections", "Journalist\nHarassment")), line=0, col="Black", col.axis="Black", col.ticks="Black", cex.axis=1.2, tick=F, las=2)
abline(v=0, lty=2)
legend("topright",bty="n",lwd=3,col=c("darkgray","black"),c("Maximum Observed Value","Minimum Observed Value"), cex=1.1)

unit=1
segments(x0=meharjrn.ish[1,2], y0=unit-adjust, x1=meharjrn.ish[1,3], y1=unit-adjust, lwd=1.5, bg="black", col="black")
points(x=meharjrn.ish[1,1], y=unit-adjust, pch=21, bg="black", col="black", cex=1.5)
segments(x0=meharjrn.ish[2,2], y0=unit+adjust, x1=meharjrn.ish[2,3], y1=unit+adjust, lwd=1.5, bg="darkgray", col="darkgray")
points(x=meharjrn.ish[2,1], y=unit+adjust, pch=21, bg="darkgray", col="darkgray", cex=1.5)

unit=2
segments(x0=frefair.ish[1,2], y0=unit-adjust, x1=frefair.ish[1,3], y1=unit-adjust, lwd=1.5, bg="black", col="black")
points(x=frefair.ish[1,1], y=unit-adjust, pch=21, bg="black", col="black", cex=1.5)
segments(x0=frefair.ish[2,2], y0=unit+adjust, x1=frefair.ish[2,3], y1=unit+adjust, lwd=1.5, bg="darkgray", col="darkgray")
points(x=frefair.ish[2,1], y=unit+adjust, pch=21, bg="darkgray", col="darkgray", cex=1.5)

unit=3
segments(x0=legcon.ish[1,2], y0=unit-adjust, x1=legcon.ish[1,3], y1=unit-adjust, lwd=1.5, bg="black", col="black")
points(x=legcon.ish[1,1], y=unit-adjust, pch=21, bg="black", col="black", cex=1.5)
segments(x0=legcon.ish[2,2], y0=unit+adjust, x1=legcon.ish[2,3], y1=unit+adjust, lwd=1.5, bg="darkgray", col="darkgray")
points(x=legcon.ish[2,1], y=unit+adjust, pch=21, bg="darkgray", col="darkgray", cex=1.5)

unit=4
segments(x0=rol.ish[1,2], y0=unit-adjust, x1=rol.ish[1,3], y1=unit-adjust, lwd=1.5, bg="black", col="black")
points(x=rol.ish[1,1], y=unit-adjust, pch=21, bg="black", col="black", cex=1.5)
segments(x0=rol.ish[2,2], y0=unit+adjust, x1=rol.ish[2,3], y1=unit+adjust, lwd=1.5, bg="darkgray", col="darkgray")
points(x=rol.ish[2,1], y=unit+adjust, pch=21, bg="darkgray", col="darkgray", cex=1.5)

unit=5
segments(x0=frassoc.ish[1,2], y0=unit-adjust, x1=frassoc.ish[1,3], y1=unit-adjust, lwd=1.5, bg="black", col="black")
points(x=frassoc.ish[1,1], y=unit-adjust, pch=21, bg="black", col="black", cex=1.5)
segments(x0=frassoc.ish[2,2], y0=unit+adjust, x1=frassoc.ish[2,3], y1=unit+adjust, lwd=1.5, bg="darkgray", col="darkgray")
points(x=frassoc.ish[2,1], y=unit+adjust, pch=21, bg="darkgray", col="darkgray", cex=1.5)

dev.print(file=paste(savedir,'backsliding_Shutdowns.pdf',sep=''), device=pdf)

# Social Media Monitoring
models=5
adjust=0.1
par(mar=c(bottom=4.0,left=6.5,top=1.0,right=0.2), lheight=0.8,xpd=F) # xpd lets plot outside graph area
plot(-1,-1, type="p", axes=F, xlab="", ylab="", xlim=c(-2,2), ylim=c(0.5,to=models+0.5), cex.lab=1.3, cex.main=1.3, main="Social Media Monitoring")
axis(1, at=seq(-2,2,by=0.5), line=0, las=0, tick=T, cex.axis=1.3)
axis(2, at=seq(1,models, by=1), label=rev(c("Freedom of\nAssociation", "Rule of Law", "Legislative\nConstraints", "Free & Fair\nElections", "Journalist\nHarassment")), line=0, col="Black", col.axis="Black", col.ticks="Black", cex.axis=1.2, tick=F, las=2)
abline(v=0, lty=2)
legend("topright",bty="n",lwd=3,col=c("darkgray","black"),c("Maximum Observed Value","Minimum Observed Value"), cex=1.1)

unit=1
segments(x0=meharjrn.smm[1,2], y0=unit-adjust, x1=meharjrn.smm[1,3], y1=unit-adjust, lwd=1.5, bg="black", col="black")
points(x=meharjrn.smm[1,1], y=unit-adjust, pch=21, bg="black", col="black", cex=1.5)
segments(x0=meharjrn.smm[2,2], y0=unit+adjust, x1=meharjrn.smm[2,3], y1=unit+adjust, lwd=1.5, bg="darkgray", col="darkgray")
points(x=meharjrn.smm[2,1], y=unit+adjust, pch=21, bg="darkgray", col="darkgray", cex=1.5)

unit=2
segments(x0=frefair.smm[1,2], y0=unit-adjust, x1=frefair.smm[1,3], y1=unit-adjust, lwd=1.5, bg="black", col="black")
points(x=frefair.smm[1,1], y=unit-adjust, pch=21, bg="black", col="black", cex=1.5)
segments(x0=frefair.smm[2,2], y0=unit+adjust, x1=frefair.smm[2,3], y1=unit+adjust, lwd=1.5, bg="darkgray", col="darkgray")
points(x=frefair.smm[2,1], y=unit+adjust, pch=21, bg="darkgray", col="darkgray", cex=1.5)

unit=3
segments(x0=legcon.smm[1,2], y0=unit-adjust, x1=legcon.smm[1,3], y1=unit-adjust, lwd=1.5, bg="black", col="black")
points(x=legcon.smm[1,1], y=unit-adjust, pch=21, bg="black", col="black", cex=1.5)
segments(x0=legcon.smm[2,2], y0=unit+adjust, x1=legcon.smm[2,3], y1=unit+adjust, lwd=1.5, bg="darkgray", col="darkgray")
points(x=legcon.smm[2,1], y=unit+adjust, pch=21, bg="darkgray", col="darkgray", cex=1.5)

unit=4
segments(x0=rol.smm[1,2], y0=unit-adjust, x1=rol.smm[1,3], y1=unit-adjust, lwd=1.5, bg="black", col="black")
points(x=rol.smm[1,1], y=unit-adjust, pch=21, bg="black", col="black", cex=1.5)
segments(x0=rol.smm[2,2], y0=unit+adjust, x1=rol.smm[2,3], y1=unit+adjust, lwd=1.5, bg="darkgray", col="darkgray")
points(x=rol.smm[2,1], y=unit+adjust, pch=21, bg="darkgray", col="darkgray", cex=1.5)

unit=5
segments(x0=frassoc.smm[1,2], y0=unit-adjust, x1=frassoc.smm[1,3], y1=unit-adjust, lwd=1.5, bg="black", col="black")
points(x=frassoc.smm[1,1], y=unit-adjust, pch=21, bg="black", col="black", cex=1.5)
segments(x0=frassoc.smm[2,2], y0=unit+adjust, x1=frassoc.smm[2,3], y1=unit+adjust, lwd=1.5, bg="darkgray", col="darkgray")
points(x=frassoc.smm[2,1], y=unit+adjust, pch=21, bg="darkgray", col="darkgray", cex=1.5)

dev.print(file=paste(savedir,'backsliding_Monitoring.pdf',sep=''), device=pdf)

# Arrests for Political Content
models=5
adjust=0.1
par(mar=c(bottom=4.0,left=6.5,top=1.0,right=0.2), lheight=0.8,xpd=F) # xpd lets plot outside graph area
plot(-1,-1, type="p", axes=F, xlab="", ylab="", xlim=c(-2,2), ylim=c(0.5,to=models+0.5), cex.lab=1.3, cex.main=1.3, main="Arrests for Content")
axis(1, at=seq(-2,2,by=0.5), line=0, las=0, tick=T, cex.axis=1.3)
axis(2, at=seq(1,models, by=1), label=rev(c("Freedom of\nAssociation", "Rule of Law", "Legislative\nConstraints", "Free & Fair\nElections", "Journalist\nHarassment")), line=0, col="Black", col.axis="Black", col.ticks="Black", cex.axis=1.2, tick=F, las=2)
abline(v=0, lty=2)
legend("topright",bty="n",lwd=3,col=c("darkgray","black"),c("Maximum Observed Value","Minimum Observed Value"), cex=1.1)

unit=1
segments(x0=meharjrn.arr[1,2], y0=unit-adjust, x1=meharjrn.arr[1,3], y1=unit-adjust, lwd=1.5, bg="black", col="black")
points(x=meharjrn.arr[1,1], y=unit-adjust, pch=21, bg="black", col="black", cex=1.5)
segments(x0=meharjrn.arr[2,2], y0=unit+adjust, x1=meharjrn.arr[2,3], y1=unit+adjust, lwd=1.5, bg="darkgray", col="darkgray")
points(x=meharjrn.arr[2,1], y=unit+adjust, pch=21, bg="darkgray", col="darkgray", cex=1.5)

unit=2
segments(x0=frefair.arr[1,2], y0=unit-adjust, x1=frefair.arr[1,3], y1=unit-adjust, lwd=1.5, bg="black", col="black")
points(x=frefair.arr[1,1], y=unit-adjust, pch=21, bg="black", col="black", cex=1.5)
segments(x0=frefair.arr[2,2], y0=unit+adjust, x1=frefair.arr[2,3], y1=unit+adjust, lwd=1.5, bg="darkgray", col="darkgray")
points(x=frefair.arr[2,1], y=unit+adjust, pch=21, bg="darkgray", col="darkgray", cex=1.5)

unit=3
segments(x0=legcon.arr[1,2], y0=unit-adjust, x1=legcon.arr[1,3], y1=unit-adjust, lwd=1.5, bg="black", col="black")
points(x=legcon.arr[1,1], y=unit-adjust, pch=21, bg="black", col="black", cex=1.5)
segments(x0=legcon.arr[2,2], y0=unit+adjust, x1=legcon.arr[2,3], y1=unit+adjust, lwd=1.5, bg="darkgray", col="darkgray")
points(x=legcon.arr[2,1], y=unit+adjust, pch=21, bg="darkgray", col="darkgray", cex=1.5)

unit=4
segments(x0=rol.arr[1,2], y0=unit-adjust, x1=rol.arr[1,3], y1=unit-adjust, lwd=1.5, bg="black", col="black")
points(x=rol.arr[1,1], y=unit-adjust, pch=21, bg="black", col="black", cex=1.5)
segments(x0=rol.arr[2,2], y0=unit+adjust, x1=rol.arr[2,3], y1=unit+adjust, lwd=1.5, bg="darkgray", col="darkgray")
points(x=rol.arr[2,1], y=unit+adjust, pch=21, bg="darkgray", col="darkgray", cex=1.5)

unit=5
segments(x0=frassoc.arr[1,2], y0=unit-adjust, x1=frassoc.arr[1,3], y1=unit-adjust, lwd=1.5, bg="black", col="black")
points(x=frassoc.arr[1,1], y=unit-adjust, pch=21, bg="black", col="black", cex=1.5)
segments(x0=frassoc.arr[2,2], y0=unit+adjust, x1=frassoc.arr[2,3], y1=unit+adjust, lwd=1.5, bg="darkgray", col="darkgray")
points(x=frassoc.arr[2,1], y=unit+adjust, pch=21, bg="darkgray", col="darkgray", cex=1.5)

dev.print(file=paste(savedir,'backsliding_Arrests.pdf',sep=''), device=pdf)
